Hydrodynamic simulations of metal ablation by femtosecond laser irradiation 
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Ablation of Cu and Al targets has been performed with 170 fs laser pulses in the intensity range 
of 10 12 -10 14 W/cm 2 . We compare the measured removal depth with ID hydrodynamic simulations. 
The electron-ion temperature decoupling is taken into account using the standard two-temperature 
model. The influence of the early heat transfer by electronic thermal conduction on hydrodynamic 
material expansion and mechanical behavior is investigated. A good agreement between experimen- 
tal and numerical matter ablation rates shows the importance of including solid-to- vapor evolution 
of the metal in the current modeling of the laser matter interaction. 

PACS numbers: 61.80.Az, 72.15.Cz, 79.20.Ap 
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I. INTRODUCTION 

Ultrafast phenomena driven by subpicosecond laser 
pulses have been the subject of thorough investigation for 
many years in order to explain the ablation of solid mate- 
rials^ 2 -^. As opposed to the laser-dielectric interaction 
where thermal and athermal ablation regimes probably 
takes place^i, the laser-metal interaction is mainly gov- 
erned by the thermal ablation one££. The laser energy is 
absorbed by the free electrons first. The pulse duration 
being shorter than the electron-phonon relaxation time 
(je-ph ~ 10 _12 s), electrons and ions subsystems must 
be considered separately. The "two-temperature model" 
(TTM) describes the thermal transport of energy by the 
electrons and the energy transfer from the electrons to 
the lattice^. 

Numerous theoretical and experimental previous works 
have been devoted to the study of the matter ablation 
with a single laser pulse. Experimental irradiation con- 
ditions in current applications are largely investigated to 
optimize the ablation rate: pulse duration^, fluenceii 
and background ga a 12 i 13 . However, a complete view in- 
cluding all the relevant physical mechanisms is still lack- 
ing. To get a better understanding of the ablation process 
and to extend the results into situations not covered by 
the experiments, two kinds of investigations are put at 
work : (i) a complete identification of the various physi- 
cal mechanisms responsible for the material removal from 
the surface, (ii) an evaluation of the impact of these var- 
ious processes on the amount of ablated matter. In the 
works previously addressed, few calculations are able to 
provide a direct comparison with experiments. Most of 
them are focused on thermal transport and a more de- 
tailed description of the physical processes occuring in 
the material has to be incorporated to really describe 
the whole ablation process. Among these, works based 
on hydrodynamic modeling^ 2 - ' 14 ' 15 have been recently as- 
sociated to the TTM to describe ablation process. To 
overcome the drawbacks of a material fluid treatment, 
a mechanical extension of the TTM has been proposed 
to model the ultrafast thermomelasticity behavior of a 



metal filmic. Works based on molecular dynamics allow 
to access to the influence of the ultrafast energy deposi- 
tion on the thermal and mechanical evolution properties 
of the materia l 17 ' 18 . With a different approach, other au- 
thors have performed a microscopic analysis of the mech- 
anisms of ultrashort pulse laser melting by means of a 
hybrid molecular-dynamic and fluid modelisation 3,4 . 

From all these investigations, it appears that none of 
these effects may be neglegted to reproduce the features 
of the damage resulting from an ultrashort laser irradi- 
ation. Moreover, there is a lack of investigations which 
combine experimental and theoretical results so that cur- 
rent models are still questionnable. In the present simula- 
tions, the TTM provides energy deposition in an expand- 
ing material intimately correlated with the processes gov- 
erning the ablation in the ultrashort pulse case, which is 
a specificity of our hydrodynamic approach. Simulation 
results give useful insight into the presented experiment 
data. 

Transport properties of electrons are not very well un- 
derstood in noncquilibrium electron-ion systems. How- 
ever, the comprehension of these phenomena in the con- 
text of ultrafast interaction is essential not only for fun- 
damental purposes but also for micromachining applica- 
tions. A precise description of the effect of the electronic 
temperature on the absorption seems to be still unset- 
tle d 19 ' 20 ' 21 , and it has not been taken into account in 
the presented calculations. Nevertheless, the model em- 
ployed in this work uses a large set of current available 
data. 

Obviously, numerical calculations are always requir- 
ing additional information such as electron-phonon or 
electron-electron relaxation times, which may be ex- 
tracted, from experimental dat a 22 ' 23 . Reciprocally, com- 
parison between simulations and experiments allows to 
validate physical data introduced into the theoretical 
models. For instance, we shall see in the following that 
the measure of the pressure variation with time inside 
the material would be very informative. These data, 
however, are difficult to measure with a high accuracy. 
By contrast, experimental measure of laser ablation rate 
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seems to be easier to be obtained and compared with nu- 
merical simulations. 

In this article, numerical and experimental results on 
ablated matter are reported. For this purpose, the TTM 
is inserted into a hydrodynamic code in order to de- 
scribe the material removal. First we detail the physical 
processes which are taken into account in our computa- 
tions within the framework of hydrodynamical modeling. 
Then, the effects caused by relaxation processes on the 
evolution of shock waves are examined. We next present 
the experiments which have been performed to obtain 
ablation depth measurements as a function of the laser 
fluence. Finally, our discussions are based on results of 
numerical simulations on Cu and Al samples compared 
with specific experimental measurements. 



II. THEORETICAL MODEL 

To represent numerically the interaction between the 
laser and the metallic target, we used a ID Lagrangian 
hydrodynamic code 24 . Solving the Hclmholtz wave equa- 
tion permitted us to determine the electromagnetic field 
through the region illuminated by the laser. The de- 
posited energy is then deduced using the Joule-Lenz law. 
A precise result for the absorption from an arbitrary 
medium can be obtained from the direct solution of the 
equation for the electromagnetic field. Let us consider 
a planar wave propagating along z axis. We write the 
following Hclmholtz equation for the complex amplitude 
of the electric field E z with frequency lj : 
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1 + i — v (T, p) j -jE z = (1) 
ui J cr 



where cr(T, p) is the complex conductivity and c is the 
light speed. The function a is calculated as a function 
of the density p(z) and the temperature T{z). The rela- 
tive permittivity of the medium is supposed to be equal 
to the unity in the case of a metal target. These simu- 
lations need accurate data such as transport coefficients 
in solids, liquids, vapors and dense or diluted plasmas, 
specially refractive indices 2 ^, electric and thermal con- 
ductivities^, and mechanical properties such as mate- 
rial strength. The evolution of the irradiated target is 
governed by the fluid hydrodynamic and heat diffusion 
equations connected with a multi-phase equation of state 
(EOS). Thermodynamic functions that realistically de- 
scribe characteristics of metal properties in various parts 
of phase diagram are needed. A such set of different 
metal EOS is generated by means of a numerical tool 
developped by the authors of the reference 2 ^. As an il- 
lustration of the EOS used, Fig. [T] displays isobars in the 
phase diagram temperature-specific energy of the copper 
for a wide range of pressure. Such data reveal the depen- 
dence on the thermodynamic properties of the melting 
and vaporisation processes. 

The mechanical propagation of shock waves and frac- 
ture are also simulated. Elastoplasticity is described 



by a strain rate independent model (relevant to high 
strain rate conditions at high pressure, typically beyond 
10 GPa) which accounts for pressure, temperature and 
strain dependent yield strength and shear modulus 28 . 
Laser induced stresses are the combination of the hy- 
drostatic pressure and the response to the shearing de- 
formation. In the temperature range of interest here, the 
effect of radiative energy transport on the hydrodynamic 
motion is negligible. Hence we ignore energy transport 
by radiation when solving the hydrodynamic equations. 
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FIG. 1: Representation of isobars in a phase diagram of the 
copper EOS in the region of phase transitions. 

Ultrashort laser irradiation and the associated ionic 
and electronic temperature decoupling require to intro- 
duce specific electronic parameters. We assume that free 
electrons form a thermal distribution during the interac- 
tion and use the Fermi-Dirac distribution to determine 
the electron properties (energy, pressure and heat capac- 
ity) as a function of the density and the temperature 2 ^. 
The evanescent electromagnetic wave is absorbed by the 
electrons. In our range of intensity, they are quickly 
heated at a temperature of few eV. During and after the 
pulse, the energy diffuses among the electrons and is then 
transferred to the lattice. Classical heat diffusion plays 
a significant role as soon as a thermal gradient occurs 
in the material. Diffusion processes are described by the 
following equations : 

p C e 



pa 



dt 
dt 



V (K e VT e ) - 7 (T e - T;) + S (2) 
V(iGVT 4 )+ 7 (r e -T 4 ) (3) 



where T, C and K are the temperature, the specific heat 
and the thermal conductivity respectively. Indices "e" 
and "i" stand for electron and ion species, p is the mass 
density. 7 characterizes the rate of energy exchange be- 
tween the two subsystems and S is the space and time 
dependent laser source term determined by the Joule- 
Lenz law. Introduction of the TTM in a hydro-code al- 
lows us to take into account the density dependence of 
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both specific heats and conductivities. C e (T e , p) is calcu- 
lated with the Fermi model using p dependent chemical 
potentials 2 —. The electron thermal conductivity K e is 
commonly expressed in the form 3 ^: 



Kr 



(<9 2 + 0.16) 5 / 4 ((9 2 + 0.44) 



(4) 



'(0 2 +O.O92) 1 /2(02 +/ 3^ ) - 

where 9 e and 9i are electron and ion temperature normal- 
ized to the Fermi temperature (9 e = T e /T F , B l — T.JT F ) 
and a, j3 are material dependent parameter s 3 ! 31 . The lin- 
ear variation of coupling term with (T e — T t ) is classic in 
TTM : we have taken 7 = 7o as 3 x 10 16 WR" 1 !^ 3 for 
copper—, and 3 x 10 17 WK _1 m~ 3 for aluminum 2 ^. It 
must be noticed that the values of K e and 7 are subject 
to considerable uncertainty in literature^. To accurately 
describe the ultrafast response, we incorporate electronic 
pressure into the set of hydrodynamic equations. The 
system of equations for electron and ion subsystems can 
be written in the Lagrangian form : 



= -P 



dV_ 

at 



(5) 



de e dV de, 

~dt = ~ Pe ~dt 7 ~di 
du 8 . ^ . dV du 

dt dm e 1 ' dt dm 

where u is the fluid velocity, m the mass and V the spe- 
cific volume. p e , e e , pi, Si are the pressure and the spe- 
cific energy of electrons and ions. Equations ([2]) to (J5j) 
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FIG. 2: Contour plots of the stress resulting from the irradi- 
ation of a copper target by a 170 fs, 10 J/cm 2 laser pulse. 

connected with a multiphase equation of state (EOS) 2 ^ 
constitute a self-consistent description of the matter evo- 
lution. However, the pressure p(T e = Tj) provided by 
this EOS is the sum of the electronic and ionic pressure 
at the equilibrium and has to be replaced by the sum 
of these two contributions out of equilibrium. The elec- 
tronic pressure is independently determined by means of 
the standard fermion gas model. As a consequence, the 
total pressure used in the above calculations is deter- 
mined as p(T e ^ T t ) = p e (T e ) - Pe(Ti) + pi{Ti). 



III. SIMULATIONS AND ANALYSIS 

To start with, the interaction of a 10 J/cm 2 , 170 fs 
FWHM gaussian pulse at 800 nm wavelength with a cop- 
per target is investigated. Fig. [5] shows the space-time 
evolution of the induced stress. The metal surface is 
heated to a maximum of 2950 K, 30 ps after the irra- 
diation. At this time, the free surface expands in a liquid 
state with a velocity of 400 ms _1 . Due to the electron 
heating, an electronic compression wave appears at the 
end of the laser pulse. The electron-ion energy exchange 
results in a significant increase in the ionic pressure, 
which propagates inside the metal. Behind the shock, 
tensile stress occurs associated with very high strain rate 
around 10 9 s -1 . 

In order to study the sensitivity of the above results 
with respect to the physical parameters, we compare in 
Fig. [3] the time variation of the computed stress at 1 /im 
depth under standard conditions (7 = 7 o and K e — K eo ) 
with those obtained with an increased coupling factor or 
electronic conductivity. 
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FIG. 3: Time evolution of the stress in copper at 1 /im depth 
resulting from a 170 fs, 10 J/cm 2 laser pulse. Standard con- 
ditions : 7 = 70 and K e = K £Q (solid line), 7 =5 70 (dashed 
line), K e = 10 K e(l (dotted line). 

In the first hundred picoseconds, the stress growth is 
directly related to the heating depth. The increased cou- 
pling factor leads to a steeper thermal gradient and a 
lower temperature at 1 pm depth compared to the other 
situations. The resulting stress is therefore lower in this 
case. The peak of the shock wave, propagating from the 
surface, reaches the 1 pm depth at 200 ps. Increasing 
the coupling factor accelerates the energy transfer from 
the electrons to the lattice and results in higher com- 
pression. Inversely, an enhanced electronic conductiv- 
ity spreads the energy spatial profile and yields reduced 
stresses. It must be noticed that in the three cases, the 
compression is followed by a high tensile stress greater 
than the characteristic tensile strength of the material. 
Nevertheless, the loading in tension is not long enough 
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(150 ps) to induce a fracture in the mediu m 32 i 33 . The 
shock duration and intensity provide a good signature 
of the balance between the electronic diffusion and the 
electron-ion coupling. Further improvements will be dis- 
cussed in the following. 

To obtain local information on the energy transfer 
induced by a femtosecond laser irradiation, a high- 
resolution time measurement of the stress reaching the 
rear side of a thin sample could be achieve d 34 i 35 . Such ex- 
perimental records could be directly compared with our 
simulations which would led us to refine the (7, K e ) phys- 
ical values accordingly. To validate the computation, we 
performed ablation measurements and compared them to 
the current simulations using standard values of (7, K e ) 
for aluminum and copper samples. 



IV. EXPERIMENT DETAILS 

Ultrashort laser pulses are generated by an ampli- 
fied all solid-state TkSaphire laser chain. Low en- 
ergy pulses are extracted from a mode-locked oscillator 
(1.6 nJ/pulse, 80 MHz, 800 nm, 120 fs). The pulses are 
then injected into an amplifying chain including : an op- 
tical pulse stretcher, a regenerative amplifier associated 
with a two-pass amplifier using a 20 W Nd:YLF laser as 
pumping source, and a pulse compressor. Linearly po- 
larized pulses with wavelength centered around 800 nm, 
an energy of 1.5 mJ at 1 kHz repetition rate and typi- 
cal duration of 170 fs were delivered. The samples are 
mounted on a three-motorized-axes system with 0.5 /zm 
accuracy. Experiments are performed in the image plane 
of an aperture placed before the objective. The result- 
ing spatial laser profile on sample is "top hat" so that 
borders and spurious conical drilling effects are reduced. 
Focusing objectives of 25 mm or 10 mm focal lengths to 
obtain fluences in a range of 0.5 to 35 J/cm 2 with the 
same beam size, 16 /xm in diameter, in the image plane. 

For ablation depth measurements, grooves are ma- 
chined by moving the sample 36 . The sample speed is 
adapted such that each point on the groove axis under- 
goes 8 consecutive irradiations at each target pass. The 
number of times the sample passes in front of the fixed 
beam can be adjusted. Fig. 0] shows a scanning electron 
microscopy (SEM) picture of the machined grooves on 
copper for 2, 4, 6, 8 and 10 passes. The groove depth 




FIG. 4: SEM picture of machined-groove profiles, from 2 to 
10 passes, on a copper sample. 



is measured using an optical profilometer with a 10 nm 
depth resolution. The ablation rates for each groove are 
deduced taking into account the sample speed, the rep- 
etition rate of the laser, the beam size and the number 
of passes. For each energy density, an averaged ablation 
rate is determined and the number of passes has been 
chosen to obtain reproducible results. From these exper- 
iments on copper and aluminum targets, we evaluate, for 
different fluences, an ablation depth averaged over a few 
tens of laser shots. The theoretical ablation depth is de- 
duced from ID numerical simulations using a criterion 
discussed hereafter. 



V. RESULTS AND DISCUSSION 

At moderate fluence, the propagation of the shock 
wave induced by the energy deposition on the lattice 
causes the surface expansion at very high speed and sig- 
nificant non-uniform strain rates. Simultaneously, the 
surface of the target is melted or vaporized as soon as the 
conditions of temperature and density required are ful- 
filled. High strain rates can turn the liquid region into an 
ensemble of droplets and ablation follows. This process 
is called the homogeneous nucleation^i. Unfortunately, 
quantitative values on the formation and ejection of liq- 
uid droplets are difficult to access because the interfacial 
solid-liquid microscopic properties of the nuclcation cen- 
ters are not accurately known. Nevertheless, in our sim- 
ulation we can consider that the liquid layer accelerated 
outside the target corresponds to the ablated matter. 
Large values of strain rates (10 9 s _1 ) indeed signal that 
droplet formation may occur and are sufficient to produce 
ablation. At higher fluence, the surface is strongly va- 
porized. The gas expands near the free surface and com- 
presses the internal matter. The vaporized part of the 
target added to the fraction of the above-defined liquid 
layer constitute the calculated ablated matter. Experi- 
mental results on ablation of copper and aluminum are 
compared with the numerical estimates in Fig. [5] and [5] 

Sharp numerical ablation thresholds are visible at 3 
and 0.5 J/cm 2 for Cu and Al targets respectively. At 
high fluence, the ablation saturates for both materials. 
This saturation occurs mainly for two reasons. Vaporiza- 
tion and subsequent gas expansion consume energy that 
does not contribute to the ablation process. Moreover, 
the liquid layer confinement increases as far as the gas 
expands and limits the liquid removal. 

As defined by Feit et alZ, a specific removal depth, 
i.e., depth removed per unit incident fluence could be a 
relevant parameter to estimate the ablation efficiency at 
a fixed pulse duration. Calculations of this quantity is 
carried out as a function of the laser fluence. Dashed 
curve presented in Fig. [5] indicate a maximum value of 
0.5 J/cm 2 in the aluminum case. The curve is smoother 
for copper in Fig. [5] and the maximum specific removal 
depth is reached at a fluence around 5 J/cm 2 . This point 
corresponds to the occurence of a critical behavior which 
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FIG. 5: Experimental and numerical (solid line) ablation 
depth as a function of the laser fluence on a copper target 
obtained with a 170 fs laser pulse. The dashed line shows 
evolution of the specific removal rate. 
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FIG. 6: Experimental and numerical (solid line) ablation 
depth as a function of the laser fluence on an aluminum tar- 
get obtained with a 170 fs laser pulse. The dashed line shows 
evolution of the specific removal rate. 



confirms a change in the hydrodynamic behavior. While 
the thickness of the liquid layer grows with the inci- 
dent energy, the specific removal depth rises. Afterwards, 
when the gas expansion starts, a growing part of the laser 
energy is transformed in kinetic energy and the specific 
removal depth drops. This suggests that an optimum 
material removal exists and refers to the situation when 
the surface vaporisation is limited. It appears that this 
quantity is relevant for material processing which can be 
optimized by operating at this optimal fluence. 

At low fluence, a noticeable discrepancy arises between 
the experimental and numerical results. The calculated 
ablated matter for a copper target is below the exper- 
imental results, while for an aluminum target, numeri- 
cal results are above. We suspect that electron trans- 
port properties should be further improved. It has been 
shown that a significant decrease in the electrical conduc- 



tivity may take place as a result of the electronic temper- 
ature increase 1 ^, which our model discards. However, the 
experimental measurements and the theoretical calcula- 
tions come to a reasonable agreement at higher fluence. 
As it has been shown by Fisher et ai2£, in the vicinity of 
a 800 nm wavelength, the interband absorption occuring 
in an aluminum target decreases with increasing temper- 
ature. The authors show that, with 50 fs laser pulses, the 
absorbtion coefficient presents a minimum near ablation 
threshold, at 5 x 10 13 W/cm 2 laser intensity. The evolu- 
tion of interband absorption with the temperature is not 
taken into account in our calculations. Consequently, we 
may overestimate the energy absorbed in this intensity 
region. 

For copper, the simulation overestimates the ablation 
threshold. This can be due to a deficit of physical pro- 
cess comprehension or to the inaccuracy of the parame- 
ters introduced in the model. No single value of 7 or K e 
can perfectly fit the sets of data shown in Fig. [5] and [6] 
As for the discussion of pressure presented above, one 
can think that a change in 7 or K e has a similar effect 
on the threshold fluence Fth- Therefore, it is interest- 
ing to investigate the fluence threshold dependence on 
7 and K e . A parametric analysis has been performed 
for a copper target. Fig. [7] displays the threshold fluence 
which has been obtained for different parameter couples 
(7, K e = K eo ) and (7 = 70, K e ). The deposited energy 
with lower thermal conductivity or higher coupling fac- 
tor can penetrate deeper into the material. As a conse- 
quence, Fth is lowered with respect to that of the refer- 
ence case (70, K eo ) and would be comparable with exper- 
imental data in the vicinity of the threshold. However, 
simulations performed in these conditions have shown a 
higher disagreement with experimental data at higher flu- 
ence due to an earlier expanded vapor. On the contrary, 
for a reduced 7 or an enhanced K e value, Fth increased 
and the ablation depth at high fluence is overestimated. 

Results obtained with one-temperature simulations ev- 
idence the importance of the TTM to reproduce the ex- 
perimental results. The good agreement obtained be- 
tween experimental data and simulations underlines the 
need of coupling the TTM model with a hydrodynamic 
code for ablation simulations in metals. Numerical re- 
sults presented in this paper give an overall description 
of processes occuring during ultrashort laser ablation ex- 
periments. 



VI. CONCLUSION 

In this paper, we have reported new results on the in- 
teraction of femtosecond laser pulses with metal targets 
at intensities up to 10 14 W/cm 2 . Numerical computa- 
tions were carried out by means of a ID hydrodynamic 
code describing the laser field absorption and the sub- 
sequent phase transitions of matter. Simulations were 
compared to original ablation experiments performed on 
aluminum and copper samples. The behavior of the ab- 
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FIG. 7: Calculation of the threshold fluence dependence on 
the coupling parameter 7 and the electronic conductivity K e 
for copper. Simulations are performed for different ratio of 
the standard value of one parameter remaining the second 
one constant. Consequently, the intersection of curves coin- 
cide with the value used in the calculation of ablation rate 
presented in Fig. [5] 

lation depth as a function of laser fluence confirms the 
importance of the use of specific two-temperature equa- 



tion of state and hydrodynamics. An optimum condition 
for laser fluence has been identified and could provide a 
precious information for efficient material processing. We 
have highlighted that ablation process is not only gov- 
erned by electronic diffusion but also by the high shock 
formation and propagation. The differences between ex- 
perimental and numerical results remain, however, more 
pronounced for low laser fluences. 

We took a great care to extend the metal properties 
to the nonequilibrium case. Nevertheless, inclusion of re- 
alistic material parameters, such as sophisticated band 
structure of metals or various scattering mechanisms, 
would allow calculations with more accuracy. Also, a 
full nonequilibrium treatment should take into account 
the conductivity dependence with both electron-electron 
and electron-phonon collisions. This work is in progress 
and implies an elaborate optical absorption model more 
suitable for ultrashort laser pulse. 

Simulations suggest that the in-depth stresses induced 
by an ultrashort laser pulse provide information of the 
matter dynamics in time, with which experimental pres- 
sure measurements could be compared. In particular, 
because it develops over temporal scales larger than the 
energy deposition one, the characteristic shape of the de- 
layed shock conveys information about the interaction 
physics and it should thus supply a promising way for 
exploring matter distortions upon picosecond time scales. 
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